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ABSTRACT 

We explore the parameter dependence of inner disk stress in black hole ac- 
cretion by contrasting the results of a number of simulations, all employing 3-d 
general relativistic MHD in a Schwarzschild spacetime. Five of these simula- 
tions were performed with the intrinsically conservative code HARM3D, which 
allows careful regulation of the disk aspect ratio, H/R; our simulations span a 
range in H/R from 0.06 to 0.17. We contrast these simulations with two previ- 
ously reported simulations in a Schwarzschild spacetime in order to investigate 
possible dependence of the inner disk stress on magnetic topology. In all cases, 
much care was devoted to technical issues: ensuring adequate resolution and 
azimuthal extent, and averaging only over those time-periods when the accre- 
tion flow is in approximate inflow equilibrium. We find that the time-averaged 
radial- dependence of fluid-frame electromagnetic stress is almost completely inde- 
pendent of both disk thickness and poloidal magnetic topology. It rises smoothly 
inward at all radii (exhibiting no feature associated with the ISCO) until just 
outside the event horizon, where the stress plummets to zero. Reynolds stress 
can also be significant near the ISCO and in the plunging region; the magnitude 
of this stress, however, depends on both disk thickness and magnetic topology. 
The two stresses combine to make the net angular momentum accreted per unit 
rest-mass 7-15% less than the angular momentum of the ISCO. 

Subject headings: accretion, accretion disks — black hole physics — MHD — radiative 
transfer 
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Introduction 



At the very beginning of accretion disk studies, their overall properties were analyzed 
by applying the constraints of energy and angular momentum conservation to the simplest 
reasonable approximation to their structure: they were assumed to b e time-steady and 



axisymmetric, and any interna 



1973 



Shakura &: Sunyaev 



1973 



vertical structure wa s integrated over (INovikov &: Thorne 



Page &: Thorne 



19741 ). The equation of energy conservation 



can be closed by counting the energy carried by photons to infinity, but no such ready 
closure exists for the angular momentum equation; angular momentum can be conserved 
for any rate of angular momentum transport through the disk provided it does not vary 
with radius. Consequently, it was necessary to guess the angular momentum flux in order 
to complete the solution. A convenient way to parametrize this guess is in terms of the net 
angular momentum accreted onto the black hole per accreted rest-mass, jnet- The choice 
made by the original papers, and still widely-used today, is to suppose that no stresses act 
on the flow from the innermost stable circular orbit (ISCO) inward to the event horizon; 
if so, Jnet = ^^(ISCO), the orbital angular momentum of a test-particle at the last stable 
orbit (here is the covariant four- velocity). 



Thornel (11974 ) 



However, this was never more than a heuristic guess. As remarked by 
although the zero-stress boundary condition is plausibly motivated by hydrodynamic 
reasoning — the inertia of matter inside the ISCO should always be much less than that 
in the stable-orbit portion of the disk outside the ISCO — it might well be invalid if 
magnetic fields are important. In fact, one of the crucial things we have learned in 
the years since the 1970s is that magnetic fields are, in fact, essential to accretion d ue 



to the presence of the magnetorotational instability (MRI) (IBalbus fc Hawlev 



On that basis , this traditional boundary condition has been questionec 



Gammie 



19991 ). and a number of numerical simulations of global disks ( iHawley fc Krolik 



Krolik 



99^ 



1999b 
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200ll: iRevnolds fc Armitagell2001 



Gammie et al 



2004 : 



Krolik et al 



Hawlev fc KrolikI I2OO2I: Machida fc Matsumotol 12003 



20051 ) have demonstrated that magnetic stresses near the 



ISCO and in the plunging region can be sizable. 

On the other hand, it has also been suggested that the magnitude of these inner-disk 
stresses may be a function of disk parameters, notably its thickness. This was the result, 
for example, of an argument based on a hydrodynamic model with constant sound speed 



(lAfshordi fc Paczyhski 



20031 ). Parameterizing the disk thickness in terms of the ratio of 



its density scale height H to radial position r. 



Reynolds fc FabianI (120081 ) found that the 



plunging region stress in a pseudo-Newtonian MH D simulation with HI R = 0.0 5 was 



rather smaller than in the analogous si mulations of 



H/R was 2-3 times largeiij. Similarly, 



5awlev fc KrolikI (12001 



Shafee et al. 



20021 ) in which 



(j2008al ) found that for a disk with 



H/R ^ 0.06-0.08 simulated with 3-d MHD in a Schwarzschild metric, the stress in the 
plunging region was significantly smaller than had been found in other simulations with 
aspect ratios a few times larger computed with different codes and somewhat different 
physical assumptions. 

It is the goal of this paper to explore how the inner disk stress depends on parameters, 
particularly disk thickness, but also magnetic geometry. To test the former dependence, 
we have performed a new series of fully general relativistic 3-d MHD simulations with 
aspect ratios H/R ~ 0.06, ~ 0.10, and ~ 0.17, all computed with the same code in the 
Schwarzschild metric and using appropriately scaled initial conditions. To explore the 
latter, we review results from previous relativistic disk simulations and, in particular, 
make detailed use of data from two previously reported simulations using a Schwarzschild 
spacetime. The disks in these two simulations had almost identical thickness, but in one 
case the initial magnetic field was a set of nested dipolar loops, while in the other the 



^We also use R to represent the radial coordinate, i.e. r = R. 
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initial field was entirely vertical. Before presenting the results of these simulations, we will 
also discuss the importance of a number of technical considerations — particularly having 
to do with spatial resolution and the establishment of inflow equilibrium — to obtaining 
meaningful results. 



Simulation Details 



The new simulations reported here were made using the code HARM3D, an intrinsically 
conservative metho d to solve the equat ions of 3- d MHD in an arbitra ry metric. This new 



code is described in 



Noble et al. 



(120091); see a 



the earlier axisymmetric version HARM, and 



so 



Gammie et al. 



Noble et al 



( I2OO3I ) for a description of 



( 120061 ) for additional details on 



the primitive variable solver. We employ the same methodology as before, with only a few 
exceptions. In the following summary of this code's techniques, we emphasize those points 
that are different from the previous description or particular to the simulations discussed in 
this paper. 

One of the principal aims of the present work is to study the influence of disk thickness 
H on the stress at the ISCO. We define it as the density moment in the coordinate frame: 



H 



-gpy/geelO - 7r/2| / 



'-9P 



(1) 



where f?^,^ is the metric, g is the determinant of the metric, and p is the rest-mass 
density. When the density profile follow s a Gaussian distribution with standard deviation 



Hg, H = y^Ha = 0.798Hg. As in 



Noble et al. 



( l2009l ). we regulate the thickness 



by cooling bound portions of the disk when the local temperature is greater than 
some target temperature T*(r). In terms of intensive quantities, bound matter satisfies 
{p + u + P) Ut > — p and gas has temperature above the target when (F — l)u/p > T*; 
here, P is the gas pressure, u is the internal energy density, is the fluid's 4-velocity, and 
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r is the adiabatic index of the equation of state: P = (F — 1) n. The relativistic enthalpy 
h = 1 + (n + P)/p. We set F = 5/3 throughout. The optically-thin cooling function is 
implemented by modifying the stress-energy conservation equation to include a sink term: 



-Cur 



The fluid-frame emissivity, £, and are designed so as to keep the density scale height 
at the desired value. The emissivity is the same as before but we slightly modified T* to 
include a neglected relativistic correction. The new target temperature is 



2 r 



Hir] 



(2) 



where Ry is the rela tivistic correction to the vertical component of gravity (lAbramowicz et al. 



1997 



Krolik 



1999a 



R.{r) = \[ll-a'{el-l)] , (3) 

Here, Ik and are the specific angular momentum (u^) and energy (ut) of circular time-like 
geodesies in the equator of a black hole with spin parameter a. For r < risco, h and 
remain at their ISCO values. 

All of the new simulations were performed in a Schwarzschild spacetime (a = 0) 
described in terms of Kerr-Schild coordinates and run for durations of 12000M to 15000M 
(in our units, G = c = 1, so that both time and distance have units M, the mass of the 
central black hole). In all cases, the initial condition was a hydrostatic torus, but we 
examined two varieties of this state: in one the radial coordinate of the pressure maximum 



r. 



Pmax 



35M and the inner edge r^^ = 20M; in the other , rp„„, = 25M and r ^r, = 15M. The 



former set of p arameters were chosen t o match those of 



match those of 



De Villiers et al 



Shafee et al. 



fl2008af ). the latter to 



(120031 ). In the remainder of this paper, we will refer to the 



^Note that Equatio n (|3l) corrects Equation ( 7.43) of iKrolikI (jl999al ). which propagated a 
typographical error in lAbramowicz et al.l (119971 ). 
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former set as "HR" and to the latter set as "LR" because the grid schemes used for the 
former were in general higher resolution than for the latter. The q parameter determines 
the angular velocity profile in the initial torus {Vt oc r~^). The choice of q along with the 
choice of pressure maximum and inner torus edge determines the characteristic vertical 
thickness of the initial torus. For LR simulations, q was set to 1.68, while for HR, q was 
determined by requiring the initial disk's thickness at r = r^^^^ to be equal to the run's 
target thickness. To study the effect of different disk thicknesses, we ran HR simulations for 
three different target temperatures, chosen to make H/R ~ 0.05, 0.08, and ~ 0.16; these 
were designated "Thin" , "Medium" , and "Thick" , respectively. We ran only LR simulations 
for the Thin and Medium cases. Finally, in all five simulations, the initial magnetic field 
consisted of dipole poloidal loops, with field lines following density contours, and with 
amplitude set such that the mean initial plasma (3 = 100. Turbulence was seeded by adding 
random perturbations to u at the 1% level. 

Boundary conditions were imposed through assignment of primitive variables in ghost 
zones; the primitive variables are p, u, (spatial velocity components), and (magnetic 
field components). The are spatial components of the four- velocity as seen by observers in 
the ZAMO frame; the 5* are the spatial components of the magnetic field as represented in 
the Maxwell field tensor, i.e., = r /yAn. Outflow boundary conditions were taken both 
at r = Tjnax and r = r^^^: all primitive variables are extrapolated at O^'^-order into the ghost 
zones, but is set to zero — and recalculated — whenever it points into the domain. In 
order to prevent numerical boundary effects from propagating outside the trapped surface, 
we chose rmin so that the numerical domain extended 5 to 25 cells inside the event horizon. 
When there was a cutout around the polar axis (see grid details below), reflective boundary 
conditions are imposed on the perpendicular vector components while all other quantities 
are extrapolated beyond the cutout with O^'^-order accuracy. If the cutout size is negligible, 
ghost zone values are set so as to make the variables continuous across the pole. In order to 
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gain a factor of four in computer resources, we simulated only a quarter of the azimuthal 
domain, employing periodic boundary conditions linking = and = 7r/2. 

Adequate resolution throughout the accretion flow and throughout the duration of the 
simulation is vital to ensure quantitative accuracy, particularly for magnetic effects. To 
explore the consequences of resolution effects, we used two different grid schemes, one for 
HR, the other for LR. Both were designed with an eye toward satisfying several conflicting 
criteria. On the one hand, it is always desirable to have the finest feasible spatial resolution. 
Toward that end, the grid scheme should provide at least several dozen cells per scale height 
on either side of the equatorial plane, and the poloidal cell aspect ratio should never be 
too large. In addition, the number of cells within a wavelength of the fastest-growing MRI 
mode should only occasionally fall below ^ 6, th e minimum resolution level at which the 



mode grows at the correct rate (jSano et al 



2004J ). On the other hand, to fit within existing 



computational resources, the total cell count cannot be too great nor the time step too 
small. 

Each run used a unique discretization tailored to its particular target thickness. 
HARM3D accomplishes nonuniform discretization through continuous coordinate 
transformations from an underlying uniform mesh. The finite volume equations are 
discretized with respect to points uniformly distributed in coordinates, x'^, which are placed 
nonuniformly in the spherical coordinate system r, 6, and (p. The center of cell Cijk is 
located at (a;i+i/2' ^i+i/2' ^fc+1/2)' '^here xf = Xq + iAx^ and Ax^ is a cell's extent in the 
direction. We choose t = x^, (p = x^, and 

Ti = e^'i , (4) 

so that Ar/r is the same everywhere. For A^i grid cells along the x^ axis and minimum and 
maximum radii Tmin and Tmax, Ax^ = 4r- log ( ^^aiai ) . 

The relationship between x^ and 9 is determined differently in the HR and LR 
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Noble et al. 



2007 



20091 ) 



simulations. In the LR group, we follow previous work (e.g., iGammie et al.l 12003 : 

but introduce a "cutout" or excised region around the polar axis: 



e{x'^) = (tt- 29c) + e sin [2 (^o + sx^)] 



(5) 



where 6c is the approximate size of the cutout, 6q and s control the nonlinearity of the 
transformation, and ^ is the amplitude of the nonlinear part. One drawback to this 
transformation is that if one tries to place a majority of the points within the first two scale 
heights from the equator, the minimum A6, which is always found at ^ = 7r/2, is so small 
that the time step is prohibitively small. 



Shafee et al 



teOOSal ). adapted here 



In the HR simulations, we follow the method of 
to include a cutout. In this method, the sinusoidal nonlinear term is replaced with a 
polynomial: 



IT 

2 



l + (l-0 (2x^ - 1) + - ^) (2x^-1)' 



(6) 



where n is a positive odd integer, 6c is the size of the excised region, and C, is still the 
relative amplitude of the nonlinear term. Note that near the equator, where most of the 
points are located, the linear term dominates and A6'(a;^) is nearly uniform. Periodicity of 
6 G [0, 27r] is ensured by making a periodic triangle function over [0,2]. We set n = 9 
whenever equation ([6]) is used for the runs presented here. 

Because gridscale dissipation scales with the ratio of cell dimension to the length 
scale on which physical quantities vary, cells that are far from cubical may have effective 
dissipation properties that are anisotropic. This might produce unphysical results because 
physical dissipation mechanisms are unlikely to have this property. For this reason, we 
strive to limit the degree of anisotropy in our cell shapes, although cells longer by factors of 
a few in the (/(-coordinate than in the others are acceptable because orbital shear tends to 
draw out features in the azimuthal direction. At the same time, it is important, especially 
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for small H/R simulations, to make the 6'-direction cell thickness small enough to put an 
adequate number of cells within a vertical scale height. Achieving these conflicting goals 
is easier with the HR grid scheme than with the LR scheme (as illustrated by the data in 
Table [H Our HR runs use grids with Ar : r/S.9 ^2:1 while achieving at least 60 cells per 
scale height. In LR simulations, although the aspect ratios are acceptable at an altitude 
~ H away from the equatorial plane, they are rather extreme close to the plane, and — of 
these runs — only ThinLR has a comparable number of cells per vertical scale height to HR 
runs (see Table [2]). 

The time step At = Aa;'' is set equal to 0.8 times the shortest cell crossing time for the 
fastest MHD characteristic from the previous time step. 

The defining parameters of the simulations are collected in Table [TJ In the first 
column we state the names of the runs. The body of the name corresponds to the disk 
thickness, the suffix refers to the grid resolution. The target aspect ratio H/R is shown 
in the second column. Note that the "radius" to which the scale height is compared is 
the radial coordinate, which is not identical to a length in any frame of reference. Since 
we regulate the temperature, not H directly, the actual value will not in general exactly 
coincide with the target. The observed value of H/R is given in Table [2] and is obtained 
by averaging over the time interval shown in its last column (whose origin is discussed in 
§ 13. 2p and radially- averaging from the ISCO (r = 6M ) to the run's r^^^^. The third column, 
cell-count, is self-explanatory. Columns 4 to 8 define the poloidal discretization. Columns 9 
and 10 state the typical cell aspect ratio at the equator and at two scale heights from the 
equator, and Column 11 lists the size of the smallest poloidal extent — A^ — of a cell, which 
always occurs at 6* = 7r/2. 

We will also analyze two older simulations, both done with the ge n eral r elativistic 



3-d MHD code GRMHD (originally described in 



De Villiers fc HawlevI tmm . Both 
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began with the gas in a hydrostatic torus. The first, called "KDOc," was described in 



Hawley fc KrolikI (j2006[ ) and had initial co nditions similar t o the LR models presented here. 



The other, called "VDO" and analyzed in 



Beckwith et al 



(j2009[ ). began with a constant 



intensity vertical magnetic field filling the annulus 35M < r < 55M and running through 
an initial torus with a pressure maximum at r = 40M. Unlike the HARM3D simulations 
there was no explicit cooling function included. Rather, the internal energy equation was 
evolved and the only heating included was shock heating captured by an artificial viscosity 
term. The resulting H/ R values are given in Table [2] and, in terms of the descriptive terms 
used here, qualify these simulations as "thick." 



3. Quality and Consistency Checks 

We are investigating the stress at the ISCO in a steady state accretion disk with a 
focus on the role of vertical scale height. To model the accretion disk system, we use the 
equations of conservative ideal MHD along with an ad hoc cooling function to control the 
disk's thickness. We can, however, run only a discrete set of simulations, evolved for a 
limited time from particular initial conditions, using a grid with restricted spatial extent 
and modest resolution. The conclusions we obtain must be assessed within the context of 
the limitations of those simulations. 

In this section we consider the effects of some of these numerical limitations by 
developing several quantitative diagnostics to measure their possible significance. We begin 
by considering the adequacy of the grid resolution, and then look at how closely the inner 
disk approximates a statistical steady state. Finally, since this study focuses on the effect 
of disk scale height for the ISCO stress, we check to see how well we are able to control that 
variable. 



-VI- 



3.1. Resolution 

Inadequate resolution can cause a number of numerical artifacts. For example, the 
growth rate of the underlying MRI can be suppre ssed if there are fewer than ^ 6 zones 



2004h . The MRI produces 



within a wavelength of the fastest growing mode ( ISano et al. 
turbulence, but only a small range in wavenumber space can be captured, possibly distorting 
the properties of that turbulence. The rate at which nonlinear mode-mode couplings 
transfer energy from large-scale motions to small may be altered. And resolution that is 
too coarse may drive an artificially large rate of magnetic numerical dissipation. 

Unfortunately, there is no a 'priori standard by which we can measure whether a 
given grid scheme either exhibits excessive magnetic reconnection or improperly evaluated 
nonlinear mode couplings. Only through a numerical convergence test, in which increasingly 
better-resolved simulations give quantitatively consistent results, can one demonstrate that 
resolution artifacts are not influencing the outcomes. As a practical matter, however, one 
carries out simulations such as these at the highest feasible resolution. Lower resolution 
simulations might provide some information, but with the computational resources at our 
disposal it has not been possible to improve sufficiently upon the resolution used in the 
production simulations to establish formal convergence. 

It is possible, however, to check whether the grid resolution satisfies certain physical 
criteria, such as having a sizable number of cells per vertical scale height and sufficient 
cells across the fastest growing MRI wavelength that the linear growth of these modes is 
correctly described. Both criteria may be met in the initial conditions, but the data must 
be examined throughout the relevant volume (the main disk body) and for the duration of 
the simulation to ensure that they continued to be satisfied. 

Data on the mean number of cells per scale height are displayed in Table |2l There we 
see that in all the HR simulations, there were at least ~ 40 cells per scale height. ThinLR 
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was almost as well-resolved as the HR simulations with respect to this measure, with 30 
cells per H. MediumLR has somewhat fewer cells per scale height, ^ 18, but even with this 
number should still be able to resolve well dynamics on the scale of a fraction of a scale 
height. 

To quantify the quality of resolution of the MRI, we evaluated the parameter 
Q = Amri/ Az, where both Amri and Az are computed in the fluid frame, 

^^^^ ^ ^^W)^'^''^' 

and 

Az = dx^'ef\ (8) 

where e'^ and are contra- and co-variant tetrad systems in the local fluid-frame, 
respectively. To assure ourselves that the simulations always satisfied this criterion, we 
created animations of Q with frames every 20M in time. In Figure [H we show sample 
stills from the simulations ThinHR, ThinLR, MediumHR, and ThickHR, each exhibiting 
a vertical slice at fixed azimuthal angle. All of these simulations exhibited better than 
adequate resolution at all times: the typical number of cells across the fastest-growing mode 
in the disk body was ^ 20. Azimuthally-averaged versions of the data shown in Figured] 
display ratios of fastest-growing wavelength to cell-size greater than 12 throughout the 
entire region within two scale heights of the equatorial plane. It is important to recognize, 
however, that even in a superbly-resolved simulation there will always be the occasional 
region in which the local poloidal field strength is small simply because this is a chaotic 
system in which fields are free to have either sign. 

The importance of Q comes from the requirement that the MRI maintain turbulence 
in the face of continual dissipation. If the MRI growth rate is artificially reduced due to 
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poor resolutionjfl and field amplification by MHD turbulence is too weak, the field strength 
can become increasingly anemic. Thus it is possible for a simulation that began with ample 
cells per MRI wavelength to suddenly find itself under-resolved. The end-result is a field 
that dies away and a cessation of accretion due entirely to inadequate resolution. 

This appears to happen in the MediumLR simulation beginning at t ^ 9000M. 
MediumLR demonstrates the utility of the resolution diagnostic. Figure [2] shows two stills, 
at times 6000M and 12000M, from the resolution animation for this simulation. Initially 
well-resolved, MediumLR becomes more poorly resolved as the field strength diminishes. 
As a result, the accretion rate also decays. 

The two GRMHD simulations can be analyzed in similar fashion. By the standards of 
the HARM3D simulations, both have relatively few cells per scale height, only ^ 12-16, not 
quite as many as in MediumLR. We also show sample snapshots of Q in Figure |3l in each 
case showing the final time-step of the simulation. KDOc was performed several years ago, 
and we did not save enough snapshots to create an animation; of all those available (every 
80M from t = 8000M to t = lOOOOM), the one shown in Figure E] appears to be the least 
well- resolved. In fact, there is an indication of a secular worsening of resolution beginning 
at t ^ 9000M. In the case of VDO, we possess more data and have confirmed that at no 
time was the resolution substantially poorer than shown here. As can be seen, in neither 



^Alth ough one generally expects po or resolution to limit the strength of the MRI tur- 



bulence, 



Fromang fc Papaloizoul ( 120071 ) found that for zero net magnetic flux unstratified 



shearing box experiments, increased resolution led to decreased turbulent stress, appearing 
to converge to zero as resolution increased. It now seems that this result arises from another 
numerical limit ation , namely the lack of any length scale besides that of a grid zone. Both 



Shi et al. 



(120091 ) and 



Davis et al. 



(120091 ) find that convergence to a nonzero stress is recovered 



when vertical gravity is added. 
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case was the resolution as good as in the best of the HARM3D simulations, although VDO 
at its end-point was better than KDOc at its. 

3.2. Inflow equilibrium 

In the end, we hope to use these evolving simulations to describe time-steady accretion 
flows. Starting from our initial condition, this state can never be achieved at all radii 
because the angular momentum removed from accreting material will be transferred to 
matter at larger radius, causing that matter to move outward. Moreover, because only a 
fraction of the initial torus mass is accreted in the duration of the simulations, the radial 
surface density profile can relax to the one associated with inflow equilibrium only within a 
short distance outside the initial inner radius of the torus (20M for HR, 15M for LR). 

Nonetheless, it is possible to identify a period of time for which an approximate state of 
inflow equilibrium really does obtain over a reasonable dynamic range in radius, subject, of 
course, to the sorts of fluctuations that occur in statistically stationary turbulent systems. 
To identify that region in both time and space, we impose several tests relying on the 
conservation of mass and angular momentum. Whenever the disk is in inflow equilibrium, 
the radial fluxes of mass and angular momentum should be constant as a function of 
radius, but because of turbulent fluctuations, they are constant only in a time-averaged 
sense. The radial range of the equilibrium region is determined by the range over which the 
time-averaged values of these quantities are nearly constant. To identify the time periods in 
which equilibrium obtains, we consider the mass interior to several specific radii M(< r; t) 
and the time-dependence of the specific angular momentum accreted onto the black hole, 
jnct- In inflow equilibrium, these should be roughly constant over time. 

The mass interior to radius r is defined as the integral of the mass density over the 
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computational volume from the horizon to radial coordinate r, or 

M{<r;t)= / dr' de dcj) p . (9) 



Thus, once inflow equilibrium is established, the mass within a given radius should stay the 
same, as the amount of mass entering from outside is matched by the mass accreted onto 
the black hole at the center. 

Figures |1H6] display the history of mass inside r = lOM, 15M, and 20M for each of 
the three disk aspect ratio categories studied. In all cases, at early times the mass in the 
disk grows as accretion from the initial torus fills in rings at smaller radii. Eventually, the 
mass-interior curves level off, signaling the approach to equilibrium with respect to this 
criterion. In the case of MediumLR, the mass in the inner disk declines at late times, a 
symptom of the diminution in the accretion rate which we attribute to the artificial decay 
of the magnetic field. 

Similar data for the two GRMHD simulations is given in Figure [71 Both simulations 
reach inflow equilibrium with respect to this criterion, after ~ 3000M in the case of KDOc, 
after ~ lOOOOM in the case of VDO. 

Another test of inflow equilibrium is provided by the history of the specific angular 
momentum accreted into the black hole, jnct- This is determined by dividing the total 
angular momentum flux by the accretion rate: 

Jnet(r,t) = ^ ' 

where the brackets represent the shell integration of the bracketed quantity: 

(X) = j d9d(f)y/^X . (11) 

The quantity jnct{r,t) will be constant in r where there is inflow equilibrium. If the 
accretion disk is truly in a statistically steady state in regard to its angular momentum flow. 
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this quantity should exhibit no trends in time, varying only modestly due to fluctuations 
intrinsic to the turbulence. Figure [H] shows Jnct(nior) as a function of time for the entire 
duration of each of the flve new simulations. They display quite different behavior. Once 
accretion begins, jnct hardly varies in ThinHR. Although there are no secular trends in 
ThinLR after t = 2000M, there are two sharp drops, at t = 4000M and 12000M, and its 
fluctuations are in general rather larger than in ThinHR. The curve of jnet(^) in MediumHR 
is more or less flat after t = 5000M, but its fluctuations are almost as large as those in 
ThinLR. By contrast, jnet rises in a succession of steps in MediumLR, showing relatively 
brief periods of near-constancy. Finally, jnet(^) rises rapidly from the beginning of ThickHR 
until about t = 8000M, holding approximately steady from then until the end of the 
simulation at 13660M. 

The time-dependence of jnet in MediumLR is instructive. Overall, especially after 
lOOOOM, jnet rises toward ^^(ISCO). As discussed in §3.H the average number of grid zones 
per most unstable wavelength dropped in this simulation as the fleld weakened. Inadequate 
resolution leads to ever-weaker magnetic fleld which, in turn, results in a rise of the speciflc 
accreted angular momentum toward the ISCO value as the magnetic stress is reduced. 

Figure [9] shows jnet for the GRMHD simulations KDOc and VDO. The zero net-flux 
simulation, KDO, shows no secular trend in this quantity after 4000M, suggesting an 
inflow equilibrium. In VDO, on the other hand, the fluctuations in jnet are much stronger, 
and there also appears to be a rising trend from the beginning of the simulation up until 
t ~ 14000M, after which the trend flattens out and the fluctuations begin to diminish. 

Combining what we have seen in the mass-interior plots with those in the jnet histories, 
we deflne the averaging periods, Atave? for these simulations as the time when both criteria 
for inflow equilibrium are met. The results of this analysis are shown in the last column of 
Table [2j In two cases (MediumHR, ThickHR), the two tests single out the same periods; in 
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one case (ThinHR), the mass-interior equilibrium period is a portion of the specific angular 
momentum equilibrium period; in four other cases (ThinLR, MediumLR, KDOc, VDO), only 
a part of the period that meets the mass-interior test is also in equilibrium according to the 
specific angular momentum test. 

With the appropriate time- averaging period chosen, we can study how the time- 
averaged accretion rate varies with radius. Of the seven simulations, in only one (ThinHR) 
is the contrast in (M) = {pu'') inside r = 20M as much as 30%; in one (MediumLR) it is 
~ 10%; in all the others (ThinLR, MediumHR, ThickHR, KDOc, VDO), it is no more than 
a few percent. 

3.3. Scale-height regulation 

Finally we examine the actual time-averaged scale heights achieved in the various 
simulations; these are shown in Figure [TDl The scale-height regulation employed in the 
HARM3D simulations is quite successful at enforcing a fixed ratio H/R (except in the 
plunging region in ThickHR), but, as shown in Table [21 the actual value obtained can 
be different from the target, ~ 20% greater in the cases of ThinHR, MediumHR, and 
MediumLR, but ~ 70% greater in the case of ThinLR. Part of this consistent offset can 
be attributed to magnetic support, and part to the fact that the temperature is typically 
slightly greater than the target. The large offset in ThinLR is a consequence of its initial 
condition, in which the gas was given a thickness almost twice as great as the target. 
Although its initial mean plasma /3 was 100, when cooling compresses the disk by a sizable 
factor, the magnetic field strengthens while the gas pressure falls. As a result, much of the 
disk mass of ThinLR outside a thin midplane layer was supported magnetically, and its 
scale height was substantially increased beyond what its gas pressure alone could support. 
This fact emphasizes the importance of choosing initial conditions relatively close to the 
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expected time-averaged state. Inside the ISCO, the aspect ratio decreases shghtly because 
the inflow time is so short that the disk cannot maintain vertical hydrostatic equilibrium. 

The time-averaged scale heights of the two GRMHD simulations are shown in Figure [TTl 
This code evolves the internal energy equation, rather than the total energy equation, and 
its simulations did not employ an explicit cooling function. Rather, entropy is conserved 
except where local shock heating is captured by an artificial viscosity term. Thus the scale 
height is not controlled directly and is determined primarily by the initial condition. In 
the end, the two simulations have very similar scale-height profiles, with H/R ~ 0.15 for 
r > lOM, but declining inward, reaching ~ 0.12 at the ISCO (r = 6M) and ~ 0.06 just 
outside the horizon. The mean H/R for both is ^ 0.14. 



With this background in mind, we can now discuss the results for time-averaged stress 
in the inner disk. We will present them in two ways: in terms of the radial profile of the 
spherical shell-integrated fluid-frame electromagnetic stress, and in terms of the angular 
momentum flux and the quantity jnct defined in the introduction. 



We begin with the radial profile of the electromagnetic stress. Both for the purpose of 
highlighting the physics and for the purpose of contrasting with the Novikov-Thorne model, 
it is best to compute it in the fluid frame: 



4. 



Results: Stress in the Inner Disk 



4.1. Fluid-frame electromagnetic stress profile 



er^ef./(47r) 



(12) 
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where is the magnetic 4- vector. Each component of the vector dx^^'' = e^'^\ dx'^ represents 
the extent of a cell's dimension as measured in the fluid element's rest frame, and e^^\ is 
the orthonormal tetrad t hat transforms v e ctors in the Boyer-Lindquist coordinate frame to 
the local fluid frame (see 



Beckwith et al. 



(l2008bl ) for explicit expressions for the tetrad). 



The vector dx^ is the Boyer-Lindquist coordinate frame version of the Kerr-Schild vector 
'^^Ks ~ [0, Ar, A6', A0] {r,6,(j)), where Ar, A6, A0 are the radial, poloidal and azimuthal 



extents of our simulation's finite volume cell located at (r, 6, 



13 



The physical significance of the electromagnetic fluid-frame stress profile is that it 
describes the rate at which angular momentum is carried outward by electromagnetic fields. 
The net angular momentum accreted by the black hole is diminished to the degree that 
stresses like these convey angular momentum outward even while inflowing matter carries 



its orbital angular momentum inward. The Novikov-Thorne model assumes that t 



le stress 



Agol fc Krolik 



fl2000h 



begins to decline outside of the ISCO, reaching zero at that point, 
showed how changing that boundary condition to account for non-zero stress at the ISCO, 
i-e., jnct < M<^(ISCO), could significantly alter the shape of the stress profile. Even in the 
disk body, well outside the ISCO, a smaller jnet can lift the time-averaged stress above the 
Novikov-Thorne curve. 

Figure fT2l shows the time-averaged W^^^{r) for each of the HARM3D simulations, 
normalized by their time-averaged accretion rate. The dotted line shows the Novikov-Thorn e 



predi ction and the dot-dash curves are examples of an Agol-Krolik profile (lAgol fc Krolik 



2000l ). Remarkably, all five of the HARM3D simulations show almost identical profiles. 



differing only in very minor ways. We concentrate on the region where inflow equilibrium 



•^Care mu s t be exerci sed to proper l y evalu ate volume integrals in the fluid frame. Both 
Krolik et al.l (120051 ) and IShafee et all ( l2008al ) correctly projected the stress tensor into the 
fluid frame, but failed to similarly project the volume element. 
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applies, here r < 20M. Outside the ISCO, but inside the domain of inflow equilibrium, 
the fluid frame stress usually (but not always) lies slightly above the Novikov-Thorne 
prediction. Although the match is not perfect, it is somewhat better described by the 
Agol-Krolik model. As the flow approaches the ISCO, where the Novikov-Thorne model 
would predict that the stress begins to fall, the electromagnetic stress rises steadily. Inside 
the ISCO, in the plunging region, the stress continues to rise inward, with a slope that is 
similar to or slightly steeper than outside the ISCO. Just outside the event horizon, the 
stress falls sharply to zero: because this stress component is nothing more than the radial 
flux of angular momentum of rotation in the disk plane, when the black hole has no angular 
momentum (i.e., does not rotate), it cannot act as a source of angular momentum and the 
stress immediately outside it must go to zero. 

The corresponding profiles for the two GRMHD simulations are shown in Figure 13. 
They are very similar to one another, and qualitatively similar to, but quantitatively 
different from, the HARM3D profiles. Like the HARM3D profiles, the radial slope of the 
stress is nearly constant in the disk outside the ISCO; unlike the HARM3D profiles, in these 
two the stress rises somewhat more steeply inside the ISCO. As a result, the peak stress in 
these two simulations is ~ 2 times greater than seen in the HARM3D cases. 

In conclusion, despite the range of temperatures and scale heights covered by these 
models, there is a great similarity between the stress profiles of all the HARM3D models as 
well as KDOc and VDO. In other words, with respect to this measure of the stress, there 
appears to be no dependence on H/R whatsoever. 

Furthermore, the similarity between KDOc and VDO suggests that the change of 
magnetic topology from closed dipolar loops on the scale of the disk thickness to net 
vertical field also makes little difference to the radial variation of accretion stress. One 
possible explanation for this insensitivity is the fact that reconnection events in the corona 
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l argely decoupled the magnetic field in the inner disk of VDO from the large-scale flux 



( iBeckwith et al 



20091 ). The remaining field in the disk then has a topology not so different 
from that in the other simulations. In §51 we review the results from previously published 
simulations to explore further the possible role of field topology on ISCO stress. 



4.2. Specific accreted angular momentum 

The value of jnet? the mean angular momentum accreted per unit rest-mass, summarizes 
the net angular momentum flow in the system. It is determined by several effects. In the 
accretion disk body, the orbital angular momentum, m^, is close to the value associated with 
a circular test-particle orbit at that radius, but can be altered by an amount ~ {H / by 
radial pressure gradients, both gas and magnetic. Stresses, both electromagnetic (Maxwell) 
and fluid (Reynolds) move angular momentum through the accretion flow; to the degree 
that they have a net divergence, they can either add or remove angular momentum from the 
fluid. In the classical Novikov-Thorne model, the fluid's angular momentum is assumed to 
match the local circular orbit angular momentum at all locations outside the ISCO, but is 
fixed at the ISCO angular momentum at all smaller radii. The stresses, whether Maxwell or 
Reynolds, are constrained to be exactly what they need to be to produce this result: finite 
in the disk body, zero in the plunging region. As a result, jnet is predicted to be exactly 
M<^(ISCO). 

As can be seen from the data listed in Table [2], this is not the case in the simulations. 
We find that, over the range of thicknesses and magnetic geometries studied, jnet ranges 
from ^ 2.93 to ~ 3.21. Thus, in all these cases, the accreted angular momentum per unit 
accreted rest-mass is 7-15% below m^(ISCO) = 3.464. Interestingly, the largest jnet by 
far was seen in KDOc, a simulation performed on a comparatively coarse grid; if it were 
discounted, the depression of jnet below would be 10-15%. 
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The separate elements contributing to this departure from the Novikov-Thorne 
prediction are shown by the curves in Figures [H] and [T51 The sohd hues correspond to the 
net specific angular momentum fiux; in infiow equilibrium this should be constant with 
radius. The dashed curves show the time-averaged specific angular momentum carried by 
the accreting matter, i.e., {phvJ'u^) / (pu^). The difference between these curves and the net 
angular momentum flux is due to the electromagnetic stress. The dot-dash curves show the 
time-averaged mass- weighted mean angular momentum at each radius, {phu^) / {p) . These 
last two quantities can be compared with the corresponding to a circular orbit, shown 
by the dotted line (which is held constant at the ISCO value inside of that radius). 

For all the models, the mass-weighted mean angular momentum generally follows the 
circular orbit value outside of the ISCO, but continues to decline inside of that radius 
rather than holding steady. In the thin disk models the offset from the circular orbit value 
is small, while in the thicker cases this offset is somewhat larger. This is precisely what one 
would expect; the gas in the disk is partially supported by the outward radial decline in 
pressure, primarily magnetic, in disks with a larger H/R. 

In every case, the curve of mean angular momentum accreted by the matter lies well 
below the curve of the local mass-weighted mean angular momentum. One way to view the 
origin of this offset stems from the fact that the material in the disk is turbulent; local fluid 
elements have angular momenta that fluctuate from time to time and from place to place. 
It should be no surprise that fluid elements with angular momentum slightly smaller than 
the value that would support a local circular orbit tend to move inward faster than those 
with larger angular momentum. In other words, the mean accreted angular momentum is 
systematically biased toward lower values by orbital mechanics that sorts the fluid elements 
according to their place in the local angular momentum distribution. 

This effect can also be identified with turbulent Reynolds stress. To see this 
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identification directly, consider the equation of angular momentum equilibrium integrated 
over spherical shells and averaged in time: 

{phu^u^) + {Ml)=3,,,,{pu'), (13) 

where is the Maxwell stress. The term oc u^u^ can be broken into two pieces, one 
reflecting the advection of the mean angular momentum (weighted by enthalpy), the other 
reflecting departures from the mean. Dividing through by the mass accretion rate, the 
previous equation becomes 

{5{pun5{hu^)) + (Mp , 
^^^-^^ ^{^) = -^-t' 

where SX = X — {X) . That is, the net rate at which angular momentum is carried inward 
per unit rest-mass accreted is the local mean angular momentum reduced by the ratio of the 
total stress. Maxwell plus turbulent Reynolds, to accretion rate. Comparing this formalism 
to Figures [H] and [151 we see that the dot-dashed curves show {phu^) / (p) ~ (hu^), while 
the dashed curves show {phu^'u^) / {—{pu^)). Their offset can then be attributed to the 
turbulent Reynolds stress normalized to the accretion rate, {{6pu^)6{hu^)) /{ — {pu^)). This 
turbulent Reynolds stress can be quantitatively significant, particularly in VDO and to a 
lesser degree in ThickHR. 

The separation between the curve of the accretion-weighted mean angular momentum 
and j'nct is the electromagnetic angular momentum flux, and in every simulation but VDO it 
clearly makes the largest contribution to the outward angular momentum flux. In all the 
other cases, the only place where the Maxwell stress does not outweigh the Reynolds stress 
is in the immediate vicinity of the event horizon. There, the matter's angular momentum 
flux becomes almost exactly the total because a non-spinning black hole has no angular 
momentum to lose. For the same reason, just outside the horizon comes to exceed 

b^b^ in magnitude, and the net electromagnetic angular momentum flux turns (weakly) 
negative. 



- 25 - 



To conclude, then, all these models show values of jnct that are reduced below the 
ISCO value, regardless of H/R. The electromagnetic stress in the fluid frame hardly varies 
at all from one simulation to the next; consequently, its contribution to the net angular 
momentum flux is likewise nearly the same in all cases. Scale height does seem to have an 
effect on the run of through the disk: u^p is reduced below the circular orbit value in 
proportion to the magnitude of radial pressure support. It also appears that the Reynolds 
stress levels may be partially controlled by the scale height in the disk and partially by the 
magnetic topology. 



Review of Previous Simulations 



In this section we summarize some results from previous simulations. Some of these, 
although not designed specifically to study the influence of scale height, can nevertheless 
provide additional information about how the inner disk stress depends on other parameters, 
notably magnetic field topology. Others are more directly comparable to our HR series of 
simulations. 



We have already discussed models K DOc flKrolik et £ 



dipole loops in the initial torus, and VDO 



vertical field piercing the initial torus. 



Beckwith et al 



Beckwith et al 



.1120051) which began with simple 



2009|), a model that began with a 



( l2008a[ ) feature three simulations all 



computed in a Kerr {a/M = 0.9) spacetime with identical fluid initial states, but differing 
initial magnetic fields: in KDPg, the field configuration was nested dipolar loops similar to 
those in KDOc; in QDP, there was a pair of quadrupolar field loops above and below the 
equator whose associated currents had opposite signs; and in TDPa, the initial condition 
held only toroidal field. 



Figure 2 of 



Beckwith et al 



(l2008al ) provides data for investigating how these different 
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initial field configurations affect inner disk stress. Panel f shows the coordinate-frame 
electromagnetic angular momentum flux (shell- and time-averaged T^{EM)) as a function 
of radius. The curves for KDPg and QDPa, the large dipolar field and quadrupolar 
simulations, are nearly identical. Both curves are generally a factor of seve ral higher than 



that o f TDPa, the one whose initial field was purely toroidal. Figure 2c of 



Beckwith et al 



( 120091 ) shows the radial run of net specific angular momentum (what we call jnet here, but 
is labeled L in that figure). Both KDPg and QDPa are slightly below (by ^ 0.05) the ISCO 
value of 2.1. TDPa, on the other hand, has an average accreted specific angular momentum 
of 2.13, very nearly the exact ISCO angular momentum, when averaged over the last 
9000M in time. The similarity between quadrupole and dipole initial fields carries over to 
simulations with a Sch warzschild hole. Another GRMHD simulation, QDO (described in 



( iBeckwith et al.ll2008b( )). also has quadrupolar loops in the initial torus. It was run for 
lO^M in time, and averaging over the last 2000M in time gives a value of jnct = 3.21 for 
the specific angular momentum accreted into the hole, similar to KDOc. As a whole, these 
results suggest that the crucial distinction may be between a field with significant poloidal 
character (either with or without net fiux) and one that is only toroidal. 

There is another case where we can clearly see the significance of a particular field 



configuration. T 



( IBeckwith et al 



(see Fig. 7 of 



le ve rtical field model VDO was evolved in both two- and three-dimensions 
20091). In the prese nt context the contrast between the two is interesting 



Beckwith et al. 



20091 ). In the axisymmetric simulation the value of jnot 
shows strong fiuctuations. Between the times of lO^M and 1.5 x lO^M (the end of this 
simulation), the mean of jnet is 2.85, but with a standard deviation of 0.43. Values as low 
as ~ 1.0 are reached at particular moments. The mean value of jnot over the same interval 
in VDO is similar, 2.89, but the standard deviation is only 0.18 and the minimum value 
reached is 2.43. Axisymmetric simulations with vertical fields typically show strong MRI 
"channel modes" characterized by extended radial fiows accompanied by radial magnetic 
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field. It is the presence of those extended radial fields through the plunging region that 
provide the strong torques. Thus, the two-dimensional simulation illustrates the basic 
principle obtaining in three-dimensions, but in exaggerated form. In this particular case, 
the presence of a net vertical field (which cannot be reconnected away within the disk) also 
prevents the antidynamo effect from dissipating the turbulence, allowing the stress in the 
plunging region to remain over the full evolution. 

3oth in nume r ical te chnique and parameters. 



Shafee et al. 



f l2008af ). The Shafee et al. 



Our HARM3D simulations are very close, 
to the general relativistic MHD simulation of 
simulation was done, like ours, in the context of a non-rotating black hole. Both employed 
intrinsically conservative Godunov algorithms differing in only minor respects. Shafee et al. 
used a temperature regulation scheme based on maintaining constant entropy rather than a 
target temperature, but — compared to our ThinHR simulation — the resulting aspect ratio 
was only slightly thicker in the mean and somewhat less constant as a function of radius. 

In their initial hydrodynamic conditions, our HR simulations were almost identical to 
those of Shafee et al., differing only in the q parameter (theirs was chosen to give an initial 
state with H/R ~ 0.1, ours had variously H/R ^ 0.05, 0.08, and 0.16). However, they did 
differ in the initial state of the magnetic field. Whereas our initial magnetic field was a set 
of nested dipole loops centered on the pressure maximum, they imposed two sets of loops, 
centered on r = 28M and r = 38M , which they then perturbed randomly with ~ 50% 
fractional amplitude. 

The spatial grid used by Shafee et al. was also very similar to the one we employed for 
ThinHR. They used a grid with 512 x 128 x 32 cells, whereas ThinHR used 912 x 160 x 64 
cells. Both radial grids were logarithmic, but their grid extended to a slightly smaller radius 
than ours (50M as opposed to 70M); their radial cell size was therefore about 1.6 times 
larger than ours. The polar-angle grid scheme for ThinHR was finer than that of Shafee 
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et al. near the equator, resulting in ~ 30% more cells per scale height. The two simulations 
had identical azimuthal resolutions because, even though we used twice as many cells, our 
azimuthal extent was also twice as great (vr/2 as opposed to vr/4). Shafee et al. reported 
that in their initial condition, there were ~ 10 cells per fastest-growing MRI wavelength, 
but do not say how that number changed during the simulation. 

Despite this very close similarity, Shafee et al. arrived at a result that was quite 
different from ours. They found that the fluid-frame electromagnetic stress followed the 
Novikov-Thorne prediction very closely all the way to r = 9M, and then maintained more 
or less that amplitude all the way to r ~ 2.5M. In contrast, as Figure \T2\ shows, in ThinHR 
(and all our other simulations), the fluid-frame electromagnetic stress is tens of percent 
above the Novikov-Thorne prediction in the disk body, and rises steeply inward inside 
r = lOM, reaching a level ^ 5 times greater than the stress found in the Shafee et al. 
calculation. Similarly, although we found ^^(ISCO) — jnct — 0.33, Shafee et al. found a 
value less than half as large, only ^ 0.14. 

The origin of this contrast is uncertain. Although we have not yet performed a 
simulation with several dipolar loops, the very minor contrast between Kerr simulations 
with dipolar and quadrupolar initial field suggests that different forms of poloidal field 
are not, by themselves, significant. However, different field structures can place different 
demands on spatial resolution; as clearly demonstrated by MediumLR (and perhaps by 
KDO), inadequate resolution can lead to substantial artificial suppression of magnetic field 
strength. The figures illustrating results from MediumLR that we have shown in this paper 
all draw on data from the period during that simulation when it remained well-resolved; at 
later times, as its resolution quality failed, its electromagnetic stresses steadily weakened 
and U0(ISCO) — jnct diminished. Because it entails more small scale structure, a pair of 
dipole loops, as in the Shafee et al. initial condition, may create a turbulent magnetic field 
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more vulnerable to reconnection than our single dipolar loop, particularly when perturbed 
by 50% and studied with a grid having larger radial cells. 

It is possible that two other considerations may also play a role in creating these 
contrasting conclusions. First, their simulation ran only to a time of lOOOOM, and they 
presented no data demonstrating how well, and over what range of radii, it reached a state 
of inflow equilibrium. Because we found that a disk of this thickness takes ^ lOOOOM to 
reach equilibrium, their shorter duration may b e problematic. Second, t he Shafee et al. 



Schnittman et al. 



(120061 ) showed that the 



simulation had an azimuthal range of only 7r/4. 
characteristic azimuthal coherence length of features in full 27r global disk simulations was 
^ 1 radian. This result suggests that an azimuthal extent of only 7r/4 might misrepresent 
the MHD turbulence. 



6. Discussion 



19911 ) that 



It is almost twenty years since the first recognition (iBalbus fc Hawleyl 
magnetic fiel ds can produce significant str ess with disks; it has been ten years since the first 



suggestions (IKrolik 



1999b 



Gammie 



19991 ) that this stress could continue within the ISCO, 
with implications for the overall efficiency and luminosity of accretion disk. However, there 
remains controversy about how large these effects may be and how they depend on disk 
parameters. Because the first general relativistic MHD accretion simul ations demonstrate d 



that these stresses are significant in disks with thickness H/R ^ 0.15 (IKrolik et al.ll2005l ) 



the discussion in the last several years has centered on whether they m ight diminish with 



decreasing disk thickness (iReynolds fc Fabian 



2008 



Shafee et al 



2008bl ). In this paper we 



have carried out a series of simulations carefully designed to isolate the effect of varying 
H/R; we find in all of our simulations that the electromagnetic fluid-frame stress increases 
steeply inward almost all the way to the event horizon. Indeed, even for a contrast of a 
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factor of 3 from the thickest to the thinnest disk, we find an almost imperceptibly small 
change in the fluid-frame electromagnetic stress profile. The most natural interpretation of 
our results is that the radial distribution of electromagnetic stress depends at most only 
weakly on disk thickness. In our simulations, these stresses diminish jnet, the net angular 
momentum per unit rest-mass accreted, by ^ 10%. This re sult supports quantitatively 



the very crude qualitative argument given in 



KrolikI ( 1l999bl ) that the Alfven speed in the 



plunging region would always become marginally relativistic, more or less independent of 
the accretion rate, so that magnetic stresses there would always be significant, but would 
always be dominated by gravitational forces. 

It is also in keeping with other lines of qualitative reasoning. One might begin by 
asking, "If stirring of MHD turbulence by the MRI leads to significa nt magnetic stresses 



i n the disk body, what might change near the ISCO?" As shown by 



Gammie fc Popham 



( 119981 ). the orbital shear in the plunging region in a Kerr spacetime differs from the 
Newtonian value by only a number of order unity. Consequently, if the magnetic pressure 
continues to be smaller than the gas pressure in the disk's equatorial plane, one would 
expect linear growth of the MRI to behave in very much the same way as in the disk 
body. If the gas pressure falls relative to the magnetic pressure so that the plasma (3 drops 
below unity, thereby quenching linear growth of the MRI, then magnetic stresses are surely 
important. However, nonlinear development of the turbulence can be expected to change 
as the infall time becomes as short or shorter than the eddy turnover time, the time for the 



energy of tur 



julent motions to move from long length scales to short. Indeed, earlier work 



( iHirose et al. 



20041 ) has shown that the magnetic field becomes both much smoother and 
somewhat less tightly- wound (i.e., the radial component grows somewhat relative to the 
toroidal component) just inside the ISCO. Such a change in structure could alter both the 
rate of field amplification and the rate of field dissipation. On the basis of simple arguments 
like these, however, it is difficult to say whether these changes should lead to a larger or 
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smaller mean field intensity and therefore stress. In the limit that both ampl ification and 



dissipation bec ome weaker, fiux-freezing results in growing magnetic stresses (IKrolik 



Gammie 



1999b 



19991). Moreover, with the principal dynamics — orbital shear and a growing radial 
velocity — all acting in the equatorial plane, there is no obvious place for a dependence on 
disk thickness. 

This physical argument is bolstered by the fact that the different simulation versions 
differ only slightly. Significantly different grid schemes, contrasting initial conditions, 
and even wholly different codes make only slight differences in the outcome. Even the 
topological contrast of substantial net vertical flux versus none at all seems to change the 
stress by only a modest amount. 

On the other hand, we have also found that fluid effects in the inner disk can also 
contribute to a diminution in jnef Pressure support of the matter in the disk is, by definition, 
proportional to (H/R)"^. Consequently, the mass-weighted mean angular momentum at any 
location in the disk is smaller than the angular momentum of a test-particle at that radius 
by a comparable amount. More significantly, Reynolds stresses can reduce it further, by 
an amount that increases both with increasing disk thickness and with net magnetic flux 
trapped on the horizon. 

A few notes of caution should be injected into t his discussion, however : Previous work 



2008al ) suggests that a 



studying accretion in a Kerr (a/M = 0.9) geometry (IBeckwith et al. 
disk magnetic field that has no net poloidal content might produce weaker electromagnetic 
stress both throughout the disk and in the ISCO region. This effect should be explored more 
thoroughly in future work. We have also shown that before interpreting simulation results 
in terms of their implications for steady-state accretion, it is important to check carefully 
both that the simulation approximates inflow equilibrium and that the simulation's spatial 
grid provides adequate resolution throughout the period studied. 
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We conclude, then, that there appears to be httle evidence for a strong dependence of 
near-ISCO electromagnetic stress on either disk thickness or the net magnetic flux. Because 
disk thickness is really a function of accretion rate, substantial near-ISCO electromagnetic 
stress should be seen in black hole accretion systems whether they are accreting at a rate 
near Eddington or far below. Its weak dependence on net magnetic flux suggests that the 
impact of electromagnetic stresses should be similarly independent of external magnetic 
boundary conditions. At the same time, we also find a supplemental reduction of the net 
accreted angular momentum due to Reynolds stresses, and this depends on both H/R 
and magnetic geometry. When the Reynolds stress is weakest, so only the near-universal 
electromagnetic stress acts, jnct is reduced below m<^(ISCO) by ^ 7-10%; when the Reynolds 
stress is strongest, the reduction is as large as 15%. 

Fully quantitative conclusions, however, await several extensions of this work: to 
rotating black holes, to disks with more complex magnetic topologies, and to disks whose 
scale heights are constant, rather than proportional to radius. In the radiation-dominated 
regime (which should apply to the inner regions of accretion disks around black holes 
whenever the accret i on rat e is near Eddington, particularly when the central mass is large: 
Shakura &: SunyaevI (119731 ) ). H is constant as a function of radius. Therefore, a fiat-topped 
disk is a more realistic model for disks as we are likely to find them in Nature. All of these 
extensions should be feasible in the near-term. 
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Fig. 1. — Sample data illustrating resolution of the MRI in simulation ThinHR (top left), 
ThinLR (top right), MediumHR (bottom left), and ThickHR (bottom right). In each case, 
the region shown is a poloidal slice at = -7r/4. Regions with deep red color are well-resolved; 
those in blue are poorly-resolved. The dashed black lines show one and two scale heights 
from the midplane. 
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Fig. 2. — Sample data illustrating resolution of the MRI in simulation MediumLR while it 
was well-resolved (left, at time t = 6000M) and after it became poorly- resolved (right, at 
time t = 12000M). In both cases, the region shown is a poloidal slice. Regions with deep 
red color are well-resolved; those in blue are poorly-resolved. The dashed black lines show 
one and two scale heights from the midplane. 
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Fig. 3. — Sample data illustrating resolution of the MRI in simulation KDOc at its end-time, 
t = lOOOOM (left) and in simulation VDO at its end-time, t = 20000M (right). In both 
cases, the region shown is a poloidal slice extending 0.3 radians, ~ 2H, from the midplane. 
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Fig. 4. — Mass (normalized to the initial total mass in the disk) inside r = lOM (thick solid 
curve), r = 15M (dotted curve), and r = 20M (dashed curve) for the thin simulations. The 
three horizontal thin solid lines show 90% of the final mass for each of these radii. (Left) 
HR. (Right) LR. 
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Fig. 5. — Mass (normalized to the initial total mass in the disk) inside r = lOM (thick solid 
curve), r = 15M (dotted curve), and r = 20M (dashed curve) for the medium simulations. 
The three horizontal thin solid lines show 90% of the final mass for each of these radii. (Left) 
HR. (Right) LR. 
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Fig. 6. — Mass (normalized to the initial total mass in the disk) inside r = lOM (thick solid 
curve), r = 15M (dotted curve), and r = 20M (dashed curve) for the thick simulation. The 
three horizontal thin solid lines show 90% of the final mass for each of these radii. 




Fig. 7. — Mass interior to r = lOM (thick solid curve), r = 15M (dotted curve), and 
r = 20M (dashed curve) for the two GRMHD simulations. The three horizontal thin solid 
lines show 90% of the final mass for each of these radii. (Left) KDOc. (Right) VDO. 
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Fig. 8. — The net specific accreted angular momentum, jj^^t, as a function of time in ThinHR 
(top left), ThinLR (top right), MediumHR (bottom left), MediumLR (bottom middle), and 
ThickHR (bottom right). Note that U(^(ISCO) = 3.464 in Schwarzschild spacetime. 
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Fig. 9. — The net specific accreted angular momentum, jnct? as a function of time in KDOc 
(left) and VDO (right). Note that ^^(ISCO) = 3.464 in Schwarzschild spacetime. 
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Fig. 10. — Time-averaged scale height as a function of radial coordinate for each of the 
HARM3D simulations. (Left) ThinHR and ThinLR, (center) MediumHR and MediumLR, 
(right) ThickHR. In each case, the heavy curves correspond to HR, the light curves to LR. 
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Fig. 11. — Time-averaged scale height as a function of radial coordinate for the two GRMHD 
simulations. The thick curve is the zero net-fiux case, KDOc; the thin curve is the non-zero 
flux case, VDO. 
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Fig. 12. — Fluid-frame electromagnetic stress, normalized by that simulation's mean ac- 
cretion rate, for each of the disk aspect ratios. (Top left) ThinHR and ThinLR, (top right) 
MediumHR and MediumLR, (bottom) ThickHR. Heavy and light curves correspond respec- 
tively to HR and LR. Dotted curves show the Novikov-Thorne model's prediction, dot-dashed 
curves the prediction of the Agol-Krolik model with additional efficiency Ae = 0.015 — a 26% 
increase relative to the Novikov-Thorne efficiency e = 0.057. 
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Fig. 13. — Fluid-frame electromagnetic stress, normalized by that simulation's mean accre- 
tion rate, for the GRMHD simulations. (Left) KDOc. (Right) VDO. Curve identifications 
are as in Fig. [121 




Fig. 14. — Accreted angular momentum per unit rest-mass for each of the three aspect 
ratios simulated with HARM3D. (Top left) ThinHR and ThinLR, (top right) MediumHR and 
MediumLR, (bottom) ThickHR. Solid curves show jnet, with the heavy curves corresponding 
to HR, the light curves to LR. Dotted curves show the angular momentum of a circular 
orbit as a function of radius; inside the ISCO, it is held constant at u^iJSCO), consistent 
with the Novikov-Thorne model. Dashed curves represent the time-averaged specific angular 
momentum carried by the accreting matter, i.e., (phu^u^) / {pu''') . Dot-dashed curves are the 
time-averaged mass- weighted mean angular momentum at each radius, i.e., {phu^)/{p). 
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Fig. 15. — Accreted angular momentum per unit rest-mass for the two GRMHD simulations. 
(Left) Simulation KDOc. (Right) Simulation VDO. Curves have the same identifications as 
in Fig. [TH 



Table 1. Simulation Parameters 



Name 


Target 
H/R 


Cell Counf* 


e Grid 




ec 


Bo s 


Cell Shapely 
{9 - V2) 


Cell Shape'^-^ 
[9 = 2H/R) 


A(?,ni„/10-3 


ThinHR 


0.05 


912 X 160 X 64 


m 


0.93 


TT- 10-15 




3.0 


: 1 


: 17 


0.48 : 1 


: 2.7 


1.4 


ThinLR 


0.05 


192 X 192 X 64 


m 


0.49 


10-15 


TT 


71 : 


1 


: 74 


2.0 : 1 : 


2.1 


0.33 


McdiumHR 


0.08 


512 X 160 X 64 


m 


0.93 


TT • 10-15 




5.8 


: 1 


: 17 


0.43 : 1 


: 1.3 


1.4 


McdiumLR 


0.08 


192 X 192 X 64 


m 


0.35 


0.083 


0.95 


5.1 : 


1 


: 5.3 


2.7 : 1 : 


2.7 


4.7 


ThickHR 


0.16 


348 X 160 X 64 


m 


0.76 


0.094 




3.0 : 


1 


: 5.3 


0.87 : 1 


: 1.4 


4.6 


KDOc 


d 


192 X 192 X 64 


CHawlev & Krolik 2006) 




0.0457r 




2.4 : 


1 


: 2.3 


2.1 : 1 : 


1.9 


11 


VDO 


d 


256 X 256 X 64 


fBeckwith et al. 2009) 




O.OItt 




2.7 : 


1 


: 2.9 


2.3 : 1 : 


2.4 


8.£^ 



'^Values of H/ R used here are the actual time-averaged quantities specified in Table [2] 
'I Unregulated scale height. 
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Table 2. Simulation Results 



Simulation Actual (iJ/i?)^ jnct A^cciis(kl < ^^/^)' Atavc/ (lO^M) 



ThinHR 


0.061 


3.13 


81 


10-15 


ThinLR 


0.085 


3.07 


60 


4.5-11.5 


MediumHR 


0.10 


3.08 


103 


5-12.5 


MediumLR 


0.091 


3.10 


35 


5-8 


ThickHR 


0.17 


2.93 


74 


8-13.66 


KDOc 


0.13^ 


3.21 


24 


4-10 


VDO 


0.14 


3.06 


32 


14-20 



'^Average of H{r,t)/R in time (over each run's At^vc) and in radius 
(weighted uniformly with respect to logr) from rigco to Tp^^^. 

^Averaging period used here was 8000M — lOOOOM instead of Atave 
because only these data were saved. 

■^Number of cells per scale height averaged with respect to logr, i.e., 
/ N'{r)d\ogr/ J dlogr, where A/'(r) is the number of cell widths within 
6 G [it /2 — H{r)/ R , 11/2 + H{r)/R] and the bounds of integration are 

^ ['''hor) ''"pmax] • 



